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Gasdynamic features of slowly rotating gas axisymmetric accretion onto a grav- 
itating center are investigated. The process of flow restructuring is studied as 



o 

H 

' the angular velocity of accreting matter approaches the Keplerian angular veloc- 

C3 ' 

ity. For spherically-symmetric accretion, the conditions are found of the existence 

X 

5_i ■ of a steady-state supersonic solution for various external boundary conditions. The 

numerical study is performed on the basis of the Lax-Friedrichs-type second order 
of accuracy high-resolution numerical scheme with the implicit approximation of 
the source term in the Euler equations. 



1 Introduction. 

Accretion onto neutron stars and black holes gives the main energy supply in galactic 
X-ray sources. Among X-ray sources with high mass companions there are long-periodic 
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pulsars, whose origin is not yet clear (Nagase 1989; Lipunov 1992). Owing to high- 
speed stellar wind, the angular momentum of the falling matter is often not sufficient 
for the accretion disk formation at the level of the Alfvenic radius where the magnetic 
pressure of the neutron star is balanced by the dynamic pressure of the falling gas. The 
angular momentum acquired in this case by the neutron star from the falling gas was 
shown (Bisnovatyi-Kogan 1991) to be such that the equilibrium rotational period of the 
X-ray pulsar might be large. It is rather difficult to make estimations for the equilibrium 
rotational period because the formation is possible of outflowing streams which carry away 
the angular momentum (Illarionov & Kompaneetz 1990; Lovelace 1995). The accretion 
picture in this case is three-dimensional and its numerical modeling is very complicated. 
It was performed only for the case of accretion onto a gravitating center, showing the 
presence of the Rayleigh-Taylor instabilities (Matsuda et al. 1989). The nature of these 
instabilities is unclear and their numerical origin can not be excluded (Steinolfson et 
al. 1994). 

To investigate the formation of long-periodic pulsars, the model is necessary for the 
accretion onto a magnetized neutron star from the stellar wind in the binary. To reduce a 
3D problem to two dimensions, we can consider either the conical accretion of nonrotating 
gas or the accretion of slowly rotating gas onto a stationary star. In the first case (Koide 
et al. 1991), the average angular momentum acquired by the neutron star is zero due to 
the symmetry of the flow, so it can not be applied to long-periodic X-ray pulsars. 

Here we use the second approach. Our aim is to consider accretion onto a magnetized 
star with the Alfvenic radius Ra ^> R* (star radius) taking into account the interaction 
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between the gas and the magnetosphere and the possible outflow formation. As the first 
step, we consider accretion with the full penetration of the gas through the magnetosphere 
substituted by the inner boundary with the radius -R* and neglect magnetogasdynamic 
effects. 

Beskin & Pidoprygora (1995) presented the approximate solution for the accretion 
of nonrotating gas onto a slowly rotating black hole in the framework of the general 
theory of relativity. In this work, we consider the accretion of rotating gas in the Newton 
approximation up to the rotation velocities close to the Keplerian velocity. 

Three different modes are usually considered as constituent parts of the astrophysical 
accretion and are fairly well investigated (Lipunov 1992; Bisnovatyi-Kogan 1989). 

1. Spherically-symmetric accretion occurs if the star velocity Voo is much smaller than 
the speed of sound in an accreting matter and the angular momentum is negli- 
gible. 

2. Cylindrical accretion is realized when > Oqo with vanishing angular momentum. 

3. Accretion disk is formed if the total angular momentum of the matter is sufficient 
for its formation. 

Real accretion is, in fact, a combination of the above-mentioned modes. 

From gasdynamic viewpoint, it is of interest to investigate the process of the transition 
from regime 1 to regime 3 for <C as the rotational velocity of the accreting matter 
approaches the Keplerian velocity. We consider the polytropic flow of ideal, perfect gas 
with the polytropic index 7 = 1.4. Calculations are performed using the second order of 
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accuracy in space high-resolution Lax-Friedrichs-type numerical scheme proposed by one 
of the authors (Barmin & Pogorelov 1995). Detailed description of the scheme is given by 
Barmin, Kulikovskiy & Pogorelov (1996). Recently, a three-dimensional Bondi-Hoyle ac- 
cretion was investigated for the accretor (star) radius varying from 10 to 0.02 Bondi radii, 
so that both subsonic and supersonic regimes of accretion could be realized (Ruffert 1994). 
In our study we are mainly interested in qualitative tracing the flow restructuring as its 
angular velocity increases. For this reason, the Euler equations are solved and the exter- 
nal flow is supposed to be supersonic. On the other hand, the size of the star is assumed 
to be sufficiently small; the inner, also supersonic, boundary is fixed at a finite distance 
from the gravitational center. The effect of the external boundary conditions is investi- 
gated on the existence and properties of a stationary supersonic accretion for the chosen 
computational region. 



2 Spherically-symmetric supersonic accretion. 

Let us consider the accretion gas flow between two spheres with the inner and outer radii 
equal to i?* and Ro respectively. If the flow at the distance i?o is assumed to be supersonic, 
then all flow parameters should be fixed at this boundary. The question is what parameter 
values can be imposed at a given distance for the existence of a stationary solution. If 
the flow is shockless and polytropic, the following three conservation equations must be 
satisfied: 

4nR 2 pU = M, (1) 
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(2) 



U\ 7 P GM 

Here p, p, U are, respectively, the density, the pressure and the radial velocity, R is 
the current distance from the gravitating center, M is the mass of the star. The values of 
the accretion rate M, the constant K, and the Bernoulli constant h t o are to be fixed at 
R = R . Introducing dimensionless variables with the units of velocity, pressure, density, 
and length equal to U , PqV 2 , po, and R*, respectively, and designating M = U /a , 
a o — (iPo/Po) 1 ^ 2 and S = GM/UqR*, we can rewrite system (l)-(3) in the dimensionless 
form using the same notations for nondimensional values of U, R, p, p, and a = ('jp/p) 1 ^ 2 
(indices "0" and "*" correspond to the outer and to the inner boundary, respectively): 



IP a 2 S_ 1 1 S 

T + 7^T"i? ~ 2 + ( 7 - 1)M 2 Rq' U 



PU = §, (5) 



p 1 



(6) 



P 1 iMl 

Being subsonic at infinity and supersonic at R = Rq, the flow becomes sonic U = 

Ub = clb at some point R = Rb > R*, where (Bisnovatyi-Kogan 1989) 

1 S 

Ub = a B = (7) 

2R B 



In the sonic point Eqs. (4)-(7) reduce to 

5-3 7 S 



= h 



to 



(8) 



4( 7 - 1) Rb 

If Rb > -Ro and the flow comes from infinity, the value h t o must satisfy the inequality 
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<h t0 < — -——S 

4(7 - 1)R 



(9) 



3 Accretion with rotation. 



In this section we consider the process of the axisymmetric rotating flow accretion onto 
a gravitating object. Analysis is performed on the basis of the numerical solution of the 
Euler gasdynamic equations. The gas is assumed to be ideal and perfect with the following 
caloric equation of state: e = p/(j — l)p, where e is the internal energy per unit mass 
and 7 = 1.4 is the polytropic index. The system of governing equations in the Cartesian 
coordinate system shown in Fig. 1 reads: 



d\J &E dG 

dt dx dz 
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Here e = p/{pj — 1) + p(u 2 + v 2 + w 2 )/2 is the total energy per unit volume. As far 
as the gravitational field is considered spherically symmetric, the polar computational 
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domain (R, 9) is chosen with the inner and outer radii equal to -R* and Ro, see Fig. 1. On 
normalizing the quantities of density, pressure, and velocity by po, cj*/?*, pou 2 Rl, where 
Po is the density at R = Ro and oj* is the gas angular velocity at R = /?* at the equator 
for the constant angular momentum distribution, the source term can be rewritten in the 
dimensionless form as 



x 



u 



(u 2 - v 2 ) + Sx 



sin 9 



2uv 



uw + Sx 



cos 9 



(11) 



(e + p)u _ u sin 9 + w cos 9 

+ Sx m 

p R 2 

Here S = GM/u 2 Rl. The form of the other vector components in Eq. (10) remains 
unchanged. 

From now on we consider only nondimensional parameters. 

The following procedure is used to construct the initial and boundary conditions. 

1. Introducing dimensionless parameter a = Uo/Uk*, where Uq is the radial velocity 
on the outer boundary and Uk* = (GM/ R*) 1 / 2 is the Keplerian velocity on the inner 
boundary, we fix the values on the outer boundary as if the flow is spherically symmetric: 

po = 1; U = aS l/2 ; Po = ^S/^M 2 , W = 0, 

where Wo is the ^-component of the velocity and M = Uo/ao- 

Corresponding parameter values inside the computational region are adopted equal to 
those at the boundary. 
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We can find the entropy and the total enthalpy dimensionless values as follows: 



^ = K; hto = S 
Po 



2 / i n \_ 

a \{ 1 -l)M 2 + 2) ifc 



2. The following initial distribution of the angular velocity is assumed: 
1 1 

uj = — = . 2 q if -R > 20 (constant angular momentum) 

uo = 0.0025 if R < 20 (constant angular velocity) 

The y-component v (normal to the plane of Fig. 1) of the velocity v is then v = uox. 

3. The values of p and p in the whole computational region are modified. Assuming 
U(R, 9) = U and W(R, 9) = 0, we find the new pressure and density values from the 
formulas 

7 p U 2 v 2 S _ p 
7 — 1 p 2 2 R pi 

After that all values on the external boundary are fixed because it is supersonic. 

No boundary conditions are necessary on the inner circle, since the flow through it is 

supersonic. 



4 Numerical scheme. 

Point clustering towards the internal boundary circle is made to obtain a sufficiently fine 
flow resolution in the vicinity of the gravitating center. The following formula is used: 

e fc _ i 

R = R* + (Ro — i?*) — — — 
with the clustering parameter (5. 
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If we introduce a polar mesh 

£, = (Z-1)A£, 1 = 1,2,. ..,L; 6 n = (n-2.5)A6, n = 1, 2, . . . , N; 
R t = R(d), A£ = 1/(L - 1), AO = n/(2N - 8) 

with the center in the accretor position, then for each cell the system (1) in the finite- 
volume formulation acquires the form 

<9Uf 

Ri ARi A9 ,n + (R l+1/2 E l+1/2n + R l _ 1 / 2 E l _ 1/2n )A9 + 

(Ef n+1/2 + E* n _ 1/2 ) ARi + R t AR t AO Hf^ 1 = (12) 

Here ARi = i^+1/2 — Ri-1/2 and E is the flux normal to the boundary, defined as: 

E = niE + n 2 G, 

where n = (711,712) is a unit outward vector normal to the cell surface. 

The assumption is made of a piecewise-parabolic distribution of the primitive gas- 
dynamic parameters q inside the cells to specify values on their boundaries and slope 
delimiters are used to obtain the nonoscillatory property. 

The averaged slope inside the cell in the radial direction is determined as (Sawada et 
al. 1989) 

A qi = [(2A J R,_i + AR^fi + (2AR l+1 + ARi)v\ARi f x , 

where 

X = AR l+1 + A Rl + AR^, „ = - * qi+1 /\ - , v - 6qi ~^ 



AR l+1 + ARi ARi + Ai?i_i 
The left and the right boundary value are then estimated as 

ft+1/2 = Qi + \ Mi, qf-1/2 = Qi-l,Mh 
9 



Aqi = mmmod(Aq h 2/j,ARi,2uARi) 



Here 



minmod(a, b, c) = < 



sign(a) min(|a|, |6|, |c|) if be > 







otherwise 



Found the values of q L and q R on the both sides of the cell boundary, the ap- 
propriate fluxes E are calculated using the modified Lax-Friedrichs formula (Barmin 
& Pogorelov 1995): 



where the matrix R is a positive diagonal matrix with the entries equal to the spectral 
radius (the maximum of eigenvalue magnitudes) of This scheme gives a drastic 
simplification of the numerical algorithm comparing to the methods based on precise 
characteristic splitting of the Jacobian matrices, while preserving nonoscillatory property. 
Having the second order of accuracy, it is less dissipative than the original Lax-Friedrichs 
scheme. Similarly, the fluxes through another pair of cell surfaces are obtained. 

As seen from Eq. (12) the source term is approximated implicitly for better stability of 
the numerical scheme. A proper linearization of this term is made to realize the numerical 
procedure of obtaining the time-converging steady-state solution: 



The promotion of the solution in time is performed with the first order of accuracy by 
the formula: 



(/ + At— )(US 1 - Uj B ) = -At[(i* +1/2 E? +1 / 2 , n /J* + ^_ 1/2 E?_ 1/2>B /i*)/Ai2« + 




H fc+i = H fe + 



(U fc+1 -U fc ). 
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(£^+1/2 + E £n-l/2) /Rl&9 + H z fc J. 

for t — k At, k — 0, 1, . . . , and At defined by the CFL condition (I is the identity matrix). 

5 Analysis of numerical results. 

All results presented in this section were obtained in the ring region with the dimensionless 
inner and outer circle radii being i?* = 1 and Rq = 100, with 56 cells in the angular 
direction and 104 cells in the radial direction, and with the clustering parameter (5 — 4. 
The calculations were performed until their full time-convergence in a quarter of the ring 
with the appropriate reflection conditions applied in the planes of the flow symmetry. 

The way of presentation the results obtained is the following. In Figs. 2-18 the isolines 
of different gasdynamic parameters and the streamlines of the flow are presented in the 
lower and upper parts of the figure divided by the rotational axis. Figures 2-3, 7-9, 
and 13-15 correspond to the whole computational region. In these figures 18 isolines 
are presented with the constant step between the maximum and the minimum value of 
functions indicated in the corners, that is, the isoline value for any function / can be found 
by the formula f\ = / min + i x (/ max — / m in)/19. In the regions with no captions in the 
corner, the streamlines are shown. Figures 4-6, 10-12, and 16-18 present the magnified 
central part (50 computational zones) of the corresponding figures related to the whole 
computational region. 

The following dimensionless parameters are chosen, at first, to study the accretion flow 
behavior for the case of slow rotation (this choice is consistent with the considerations 
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from Section 2): 

a = 0.1; 7 = 1.4; M = 2; 5 = 250 

The steady state results for this case are presented in Figs. 2-6 (the axis of rotation is 
aligned with the ,2-axis (see Fig. 1). In Fig. 2, the isolines of the pressure and the density 
decimal logarithm are presented above and below the rotational axis, respectively, in the 
whole computational region. The isolines of the velocity components v and W (^-direction) 
are not shown in the full computational region since their main variation is located in the 
vicinity of the inner boundary. The isolines of the velocity component U (radial direction) 
and the streamlines are given in Fig. 3. 

In Figs. 4-6 the same lines are shown in the magnified subregion close to the accreting 
center (50 computational zones). It is clearly seen in these figures that for long distances 
from the accreting center and far enough from the rotational axis slow rotation does not 
affect the flow significantly and it is, in fact, the superposition of a spherically symmetric 
accretion and an axisymmetric constant angular momentum rotation. Closer to the z- 
axis and to the accreting center, however, the pressure and the density are greater at the 
equator than near the poles. The deviation from the purely constant angular momentum 
rotation appears only in the vicinity of the rotational axis and near the accreting object. 
The streamlines and the radial velocity U isolines behave quite like in the spherically 
symmetric case. 

The similar pictures of the flow are presented, for S = 25 with the rest of dimensionless 
parameters unchanged, in Figs. 7-12. The effect of rotation in this case is quite definite 
in the whole region surrounding the rotational axis. The gas displacement from the poles 
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is much more pronounced comparing with the variant of S — 250 even at large distances 
from the accretion center, see Fig. 7. The streamlines deflect from the poles to the 
equator region (Fig. 12) and the size of the domain with substantial values of W is larger 
(Figs. 8, 11). 

The picture of the flow for the rotational velocity close to the Keplerian velocity 
(S = 1.5) is shown in Figs. 13-18. Almost all the gas is removed in this case from the 
pole regions. The density near the equator is ~ 10 4 times greater than that near the poles 
and almost all accreting matter is involved in the motion towards 6 = ir/2. The structure 
of the flow is very close to the picture of the accretion disk formation, see Fig. 18 showing 
the streamline distribution. The regions with the dense distribution of the pressure and 
density isolines (Fig. 16) represent the oblique transverse shock where supersonic flow of 
the gas in ^-direction, declined by the centrifugal force from the poles to the equator, 
decelerates to become zero at 9 = n/2. Behind these shocks the region of a fast rotation 
is observed (Figs. 14, 17). This region of a highly rotating dense matter will form an 
accretion disk if S < 1. 

6 Discussion. 

Numerical modeling is performed for axisymmetric rotating gas accretion onto the star. 
The second order of accuracy in space high-resolution Lax-Friedrichs-type scheme showed 
good performance for the calculation of flows with an extremely great variation of the 
pressure and density values throughout the computational region. The flow restructuring 
is investigated as the rotation velocity approaches the Keplerian velocity. The variants 
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are presented for different cases begining from a slow rotation and up to the case, for 
which the steady state centrifugal and gravitational forces near the stellar equator are 
rather close. The values of the angular momentum, however, were sufficiently small, so 
that at the inner boundary the centrifugal force was smaller than gravitational, and the 
stationary accretion picture could be attained. If the angular momentum of the falling 
matter reaches the value, for which the centrifugal force balances the gravitation before 
the matter falls onto the star (the inner boundary in our calculations), the matter stops 
near the equatorial plane and forms a disk. In the absence of viscosity the mass of this 
disk increases in time. In the presence of viscosity, which is especially efficient in the 
turbulent case, the matter in the disk looses its angular momentum and moves slowly 
towards the center forming a stationary accretion disk (Lynden-Bell 1969; Shakura 1972). 
We are interested here in the problem of the long-periodic pulsar formation, for which 
the angular momentum is not sufficient for the accretion disk origin near the Alfvenic 
surface where the dynamic pressure of the flow is balanced by the magnetic surface at the 
magnetosphere of the neutron star. To find the angular momentum of the matter falling 
onto the star during accretion from the stellar wind penetrating through the Alfvenic 
surface, we need to include into consideration the magnetic field. This work is now in 
progress. 
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Figure captions. 

Fig. 1. Computational region. 

Fig. 2. Logarithm pressure and density isolines. Full computational region, S = 250. 

Fig. 3. Streamlines and U isolines. Full computational region, S = 250. 

Fig. 4. Logarithm pressure and density isolines. Inner subregion, S = 250. 

Fig. 5. Isolines of v and W. Inner subregion, S = 250. 

Fig. 6. Streamlines and U isolines. Inner subregion, S = 250. 

Fig. 7. Logarithm pressure and density isolines. Full computational region, S = 25. 

Fig. 8. Isolines of v and W. Full computational region, S = 25. 

Fig. 9. Streamlines and U isolines. Full computational region, S = 25. 

Fig. 10. Logarithm pressure and density isolines. Inner subregion, S = 25. 

Fig. 11. Isolines of v and W. Inner subregion, S = 25. 

Fig. 12. Streamlines and U isolines. Inner subregion, S = 25. 

Fig. 13. Logarithm pressure and density isolines. Full computational region, S = 1.5. 

Fig. 14. Isolines of v and W. Full computational region, S — 1.5. 

Fig. 15. Streamlines and U isolines. Full computational region, S = 1.5. 

Fig. 16. Logarithm pressure and density isolines. Inner subregion, S — 1.5. 

Fig. 17. Isolines of v and W. Inner subregion, S = 1.5. 

Fig. 18. Streamlines and U isolines. Inner subregion, S — 1.5. 
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